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A  PRACTICAL  CONFUTATION AL  METHOD  FOR  REDUCING  A  DYNAMICAL  SYSTEM  WITH 
CONSTRAINTS  TO  AN  EQUIVALENT  SYSTEM  WITH  INDEPENDENT  COORDINATES 


William  C.  Walton,  Jr.,  and  Earl  C.  Steeves* 
NASA  Langley  Research  Center 
Langley  Station,  Hampton,  Va. 


A  new  method  is  presented  by  which  equations  of  motion  of  a  linear 
mechanical  system  can  be  derived  in  terms  of  independent  coordinates 
when  the  system  is  described  in  terms  of  coordinates  which  are  not 
independent  but  instead  are  governed  by  linear  homogeneous  equations 
of  constraint.  There  is  a  discussion  of  the  origin  in  practical  vibrations 
analysis  of  dynamical  systems  involving  equations  of  constraint. 
Methods  previously  used  for  handling  such  systems  are  discussed  and 
the  new  method  is  demonstrated  to  have  the  following  advantages: 
(1)  For  the  most  general  constraint  equations,  solution  of  the  equations 
is  reduced  in  substance  to  computing  the  eigenvalues  and  eigenvectors 
of  a  symmetric  matrix;  and  (2)  the  method  is  applicable  when  there 
are  redundancies  in  the  equations  of  constraint. 
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SECTION  I 
INTRODUCTION 


The  purpose  of  this  paper  is  to  present  a  method  by  which  equations  of  motion  of  a  linear 
mechanical  system  can  be  derived  in  terms  of  independent  coordinates  when  basic  information 
about  the  system  is  available  in  terms  of  coordinates  which  are  not  independent  but  instead 
are  governed  by  linear  homogeneous  equations  of  constraint.  Necessity  for  this  derivation 
arises  frequently  in  practical  vibration  analysis.  The  method  is  believed  to  be  new,  and  ex¬ 
perience  in  analyzing  the  vibrations  of  shells  has  convinced  the  authors  that  it  very  often 
offers  decided  advantages  over  methods  previously  used. 

There  is  a  discussion  of  the  reason  dynamical  systems  involving  constraint  equations 
arise  in  practical  analysis  of  oscillations  of  mechanical  systems.  There  follows  a  description 
of  methods  previously  used  in  dealing  with  such  systems.  Then  a  theorem  designated  the 
“zero  eigenvalues  theorem”  which  is  basic  to  the  method  of  this  paper  is  proved,  and  the 
method  is  presented.  Next,  the  result  of  the  method  of  this  paper  is  shown  to  include  the 
result  of  the  main  older  method  as  a  special  case.  Two  examples  of  application  of  the  new 
method  are  presented  and  the  paper  closes  with  a  discussion  of  numerical  considerations 
involved  in  practical  computing  with  the  method. 
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SECTION  n 
BACKGROUND 


In  conventional  analyses  of  small  forced  oscillations  of  mechanical  systems ,  the  physical 
system  is  idealized  so  that  its  configuration  at  any  instant  is  determined  by  specification  of 

a  finite  number  of  independent  coordinates  q1 ,  q2 , .  .  q^ . qN>  Then,  with  approximations 

allowable  because  of  the  assumed  smallness  of  the  oscillations,  the  Lagrangian  of  the  system 
may  be  expressed  in  the  form 

L  --  q'Mq  -  y  q'Kq  (I) 

where 

(1)  q  is  a  column  matrix  the  elements  of  which  are  the  coordinates  q  . 

n 

(2)  A  prime  denotes  the  transpose  of  a  matrix. 

(3)  M  and  K  are  constant  symmetric  matrices  of  order  N  with  M  positive  definite. 

(4)  A  dot  denotes  differentiation  with  respect  to  time. 

When  the  Lagrangian  has  the  form  shown  by  Equation  1,  application  of  Hamilton’s 
principle  gives  equations  of  motion  of  the  system  which  have  the  form 

Mq  +  Kq  =  Q  (2> 

In  Equation  2  Q  is  a  column  matrix  with  N  elements.  The  elements  of  Q  are  usually  called 
generalized  forces;  and  for  forced  oscillations,  the  problem  under  consideration,  they  are 
functions  of  time  alone.  The  generalized  forces  are  determined  by  the  following  requirement: 
Let  8  q  be  an  arbitrary  infinitesimal  variation  of  the  coordinates  composing  the  Matrix  q  . 
Then  the  work  W  done  by  the  instantaneous  forces  applied  to  the  system  when  these  forces 
are  moved  through  the  displacements  produced  by  the  variation  shall  be  given  by  the  equation 

W  =  Q'  8  q  <  3) 

Once  the  equations  of  motion  have  been  brought  to  the  form  indicated  by  Equation  2  there 
is  a  well-established  and  quite  effective  body  of  mathematical  theory  and  computational 
technique  for  determining  the  behavior  of  the  system. 

Often,  however,  it  is  much  easier  to  express  L  and  W  if  the  system  is  described  by 
coordinates  which  are  not  independent  but  which  are  governed  by  linear  homogeneous  equations 
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of  constraint.  Letting  z^,  z^,  . 
constraint  equations  take  the  form 


z 

P 


represent  such  a  set  of  coordinates,  the 


C  z  =  0 


(4) 


where  z  is  a  column  matrix  the  elements  of  which  are  the  coordinates  z  ,  and  where  C  is  a 

P 

constant  matrix  which  has  P  columns  and  is,  in  general,  rectangular. 


In  terms  of  the  coordinates  the  Lagrangian  will  generally  take  the  form 

L  =  -y  z'Mi-i-  z'Kz 


(5) 


where  M  and  K  are  symmetric  matrices  of  order  P.  The  work  W  can  be  found  in  the  form 


W  =  Z  8z 


(6) 


where  S  z  is  an  arbitrary  variation  of  z  compatible  with  the  Equations  of  Constraint  4, 
and  Z  is  a  column  matrix  with  P  elements  which  are  functions  of  time  alone.  It  is  sometimes 
possible  to  choose  the  coordinates  z^  so  that  the  Matrix  M  is  positive  semidefinite  rather 
than  positive  definite;  but  in  this  paper  such  choices  are  excluded  and  M  is  assumed  to  be 
positive  definite. 


From  what  has  been  said  it  can  be  seen  that  it  is  useful  to  know  a  systematic  procedure 
by  which  equations  of  motion  in  terms  of  independent  coordinates,  as  in  Equation  2,  can  be 
derived  starting  with  the  Lagrangian  L  and  the  work  W  in  terms  of  coordinates  governed  by 
homogeneous  equations  of  constraint  as  in  Equations  5  and  6.  It  is  the  object  of  this  paper  to 
set  forth  such  a  procedure,  but  before  doing  so  it  is  appropriate  to  discuss  briefly  how  the 
problem  has  been  dealt  with  in  the  past. 
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SECTION  m 

OLDER  METHODS 


In  the  past,  the  equations  in  terms  of  independent  coordinates  have  been  arrived  at  in  two 
ways: 

(1)  Through  consideration  of  particular  physical  or  geometrical  aspects  of  a  problem 
the  dependent  coordinates  z  are  chosen  so  as  to  impart  a  very  simple  form  to  the  equations 
of  constraint,  rendering  easy  and  obvious  the  determination  of  independent  coordinates. 

(2)  Certain  of  the  coordinates  are  selected  to  be  independent  coordinates,  and  the 
equations  of  constraint  are  then  solved  as  simultaneous  equations  to  express  the  remaining 
dependent  coordinates  in  terms  of  those  which  have  been  selected  to  be  independent. 

Under  the  first  approach  come,  for  example,  those  finite  element  methods  of  structural 
analysis  in  which  the  coordinates  of  a  free-body  element  are  displacements  and  rotations  at 
nodes.  In  such  analyses  the  equations  of  constraint  come  down  to  equalities  among  appropriate 
displacements  and  rotations  at  moving  nodes  and  equations  in  which  appropriate  displacements 
and  rotations  are  set  equal  to  zero  at  nodes  where  there  are  supposed  to  be  rigid  constraints. 
A  set  of  independent  coordinates  is  arrived  at  by  the  simple  expedient  of  using  a  single  symbol 
for  all  displacements  and  rotations  which  are  equated  at  a  node.  This  is  the  basis  of  the  now 
widely  used  procedure  of  superimposing  stiffness  matrices  or  mass  matrices  of  structural 
elements  to  arrive  at  a  stiffness  matrix  or  a  mass  matrix  of  an  entire  structure  composed 
of  the  elements  connected  together. 

As  a  basis  for  describing  the  advantages  of  the  method  to  be  presented,  the  second  type 
of  approach  will  be  discussed  formally.  In  this  approach  it  is  assumed  (usually  tacitly)  that 
the  rank  R  of  the  Matrix  C  is  equal  to  the  number  of  rows  in  C  and  that  therefore  Equation  4 
may  be  written  as 

AilQ,  +  B:(b)  =0  (7) 

where: 

(1)  A  is  an  R  by  R  nonsingular  matrix  the  columns  of  which  are  R  distinct  colums  of  C  . 

(2)  B  is  an  R  by  (P-R)  matrix  the  columns  of  which  are  those  columns  of  C  not  included 
in  A  . 


311 


AFFDL-TR-68-150 


(3)  z  ^  and  z  are  column  matrices  the  elements  of  which  are  elements  of  z 
corresponding  to  the  columns  in  A  and  B,  respectively,  and  appropriately  ordered. 


.(b) 


By  renumbering  the  coordinates  and  making  a  corresponding  rearrangement  of  the 
columns  of  C  it  can  be  contrived  that  the  first  R  columns  of  C  constitute  the  Matrix  A  and 

(cl) 

the  last  P-R  columns  of  C  constitute  the  Matrix  B.  Correspondingly,  the  elements  of  z' 
would  be  the  first  R  elements  of  z  and  the  elements  of  z  ^  the  last  P-R  elements  of  z .  For 
convenience  in  the  ensuing  discussion  it  is  assumed  that  such  a  rearrangement  has  been 
made.  However,  as  a  practical  matter,  it  is  very  important  to  note  that  actually  to  carry  out 
a  suitable  rearrangement  one  must  be  able  to  recognize  R  linearly  independent  columns  of  C  . 
This  may  not  be  easy. 


Since  A  is  nonsingular,  an  inverse  of  A  exists  and  is  unique.  Equation  7  is  satisfied 
therefore  if  —  and  only  if 


,  (  a) 


-A'1  Bz 


(b) 


(8) 


where  A-1  is  the  inverse  of  A  .  It  follows  that  the  Equations  of  Constraint  4  are  satisfied 
if  —  and  only  if 


2 


(9) 


where 


(3- 


-a"1  B 


(10) 


In  Equation  10  I  is  an  identity  matrix  of  order  P-R. 


Substitution  of  Equation  9  into  Equations  5  and  6  gives  expressions  for  the  Lagrangian  L 
and  the  work  W  in  terms  of  independent  coordinates  and  in  the  forms  shown  by  Equations  1  and 
3,  respectively.  The  ingredients  of  the  resulting  equations  are 


q  =  z(b) 

(ID 

m  =  0fi*  0 

(12) 

k  =  0'^0 

(13) 

Q-.  0'  z 

The  matrices  M  and  K  thus  derived  are  symmetric  and, 

clearly  linearly  independent,  the  Matrix  M  is  positive  definite. 

(14) 

since  the  columns  of  0  are 
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The  preceding  is  a  fair  description  of  the  textbook  method  for  handling  equations  of 
constraint.  For  emphasis,  it  is  noted  once  more  that  the  method  requires  that  the  Rank  R  of 
the  Matrix  C  be  equal  to  the  number  of  rows  of  C  and  that  one  be  able  to  recognize  R 
linearly  independent  columns  of  the  Matrix  C . 


SECTION  IV 

THE  ZERO  EIGENVALUES  THEOREM 


The  object  here  is  to  prove  a  theorem  which  is  the  foundation  of  the  method  of  this  paper. 


Consider  the  equation 

Cz  :  0  (15) 

where  C  is  a  matrix  with  any  number  of  columns  and  any  number  of  rows.  Let  P  be  the 
number  of  columns. 


Let  a  Matrix  E  be  defined  by  the  equation 

E  =  C'  C  (16) 

and  note  that  E  is  symmetric  of  order  P. 


Since  E  is  symmetric  one  can  always  find  a  square  Matrix  U  of  order  P  such  that 

U'U  =  I  (17) 


and 

U'  EU  =  X  (18) 

where  I  is  an  identity  matrix,  and  X  is  a  real  diagonal  matrix,  each  of  order  P.  Any  matrix 
U  having  these  properties  is  called  a  modal  matrix  of  E .  The  numbers  occupying  the  main 
diagonal  of  X  are  called  the  eigenvalues  of  E.  Let  Xp  represent  the  eigenvalue  at  the  inter¬ 
section  of  the  p^1  row  and  p^1  column  of  X .  Then  we  say  that  the  p^  column  of  the  modal 
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matrix  U  is  an  eigenvector  of  the  Matrix  E  corresponding  to  the  eigenvalue  X  p.  From 
Equation  17  then  the  P  columns  of  a  modal  matrix  constitute  an  orthonormal  set  of  eigen¬ 
vectors  of  E  .  For  convenience  it  is  assumed  that 

\>X>X>  ->Xn  (19) 

12  3  P 

since  the  positions  of  the  eigenvalues  can  be  reordered  simply  by  reordering  the  columns  of 

U. 


Now,  let  a  Matrix  D  be  defined  by  the  equation 

D  =  CU  (20) 

Then  Equation  18  may  be  written 

d'd  =  X  (21) 

It  is  clear  from  the  form  of  E  shown  in  the  defining  Equation  16  that  E  is  positive 
semidefinite.  It  follows  by  well-known  theorems  that  the  eigenvalues  are  positive  or  zero. 
Let  S  represent  the  number  of  them  which  are  positive.  Then  the  last  P-S  are  zero.  It  follows 
from  Equation  21  that  the  last  P-S  columns  of  D  are  null  and  the  first  S  columns  are  linearly 
independent  {in  fact,  orthogonal  one  to  another). 


Making  use  of  Equation  17,  Equation  15  may  be  written  as 

Cliu'z  s  0  (22) 

or 

Dz  =  0  (23) 

where 

z*  =  u'  z  (24) 

Further,  substitution  of  Equations  20  and  24  into  Equation  23  and  subsequent  substitution 
of  Equation  17  returns  Equation  15  uniquely. 


Therefore,  for  a  column  z  to  satisfy  Equation  15  it  is  both  necessary  and  sufficient  that 
the  column  z  defined  by  Equation  24  should  satisfy  Equation  23. 

But  T  satisfies  Equation  23  if  —  and  only  if  —  the  first  S  elements  of  z  are  zero. 
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From  Equations  17  and  24 

z  =  Uz  (25) 

Thus  we  have  the  basic  theorem  that  z  is  a  solution  of  Equation  15  if  —  and  only  if  —  z 
may  be  expressed  in  the  form 

Z  =  T  tf  (26) 

where  q  is  a  c olumn  matrix  with  P-S  arbitrary  elements  and  where  T  is  a  matrix  the  columns 
of  which  are  those  eigenvectors  in  U  corresponding  to  eigenvalues  with  the  value  zero. 


SECTION  V 

PROPOSED  METHOD 

With  the  basic  theorem  from  the  preceding  section  in  hand,  the  following  procedure  may 
be  proposed. 

GIVEN: 

(1)  K  and  M  each  constant  symmetric  matrices  of  order  P  with  M  postive  definite. 

(2)  C  ,  a  constant  matrix  with  P  columns  and  any  number  of  rows. 

(3)  Z  ,  a  column  matrix  with  P  elements  each  of  which  may  be  a  function  of  time. 

OBJECT: 

(1)  To  compute  a  Matrix  T  such  that: 

(a)  The  transformation  z  =  Tq  relates  the  dependent  coordinates  z  appearing  in 
Equations  5  and  6  to  a  set  of  independent  coordinates  q  suitable  for  use  in  Equation  2. 

(b)  The  transformation  Q  =T'Z  produces  a  Matrix  Q  suitable  for  use  in  Equation  2. 

(2)  To  compute  Matrices  K  and  M  suitable  for  use  in  Equation  2. 
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PROCEDURE: 

(1)  Compute  E  where  E  =  c'C.  E  will  be  square-symmetric  of  order  P  and  positive 
semidefinite. 

(2)  Compute  a  modal  Matrix  U  and  the  eigenvalues  X  of  the  Matrix  E .  (p  =  1,2,3,. . .  ,P.) 

r 

This  is  a  standard  operation  at  modern  computing  installations  and,  in  fact,  is  one  of  the 
most  successful  applications  of  digital  computers. 

(3)  Identify  the  columns  of  U  which  correspond  to  zero  eigenvalues.  This  step  requires 
attention  because  in  principle  one  can  fairly  question  the  possibility  of  a  rigorous  distinction 
between  finite  eigenvalues  and  eigenvalues  having  the  value  zero  when,  as  is  normal,  there  is 
roundoff  error  in  the  process  by  which  the  eigenvalues  are  computed.  The  point  is  discussed 
in  the  section  on  numerics. 

(4)  Assemble  a  matrix  the  columns  of  which  are  the  columns  of  U  corresponding  to  the 
eigenvalues  having  the  value  zero.  This  matrix  is  the  required  transformation  matrix  T  .  Its 
dimensions  are  P  by  (P-S)  where  S  is  the  number  of  finite  eigenvalues  of  E . 

(5)  Compute  M  and  K  by  the  formulas  M  =T(  MT  and  K  =  Tx  KT  .  M  and  K  will 
each  be  symmetric  and  M  will  be  positive  definite,  provided  sufficient  numerical  significance 
has  been  carried  in  all  computations. 


316 


AFFDL-TR-68-150 


SECTION  VI 

GENERALIZATION  OF  THE  TRANSFORMATION 


Equation  26  gives  the  most  general  solution  to  the  Equations  of  Constraint  4.  Since  the 
Matrix  q  appearing  in  Equation  26  is  completely  arbitrary,  the  solution  can  just  as  well  be 
stated  in  the  form 

z  -  TH  q*  (27) 

where  H  is  any  nonsingular  square  matrix  of  order  (P-S). 


Therefore  letting 

T  -  TH  (28) 

the  Matrix  T  may  be  used  as  the  transformation  matrix  in  place  of  T . 

In  order  that  a  Matrix  T  may  be  written  as  in  Equation  28,  for  some  Matrix  H  .  It  is  both 
necessary  and  sufficient  that  the  columns  of  T  constitute  a  set  of  linearly  independent  eigen¬ 
vectors  of  E  corresponding  to  the  eigenvalues  of  E  which  have  the  value  zero.  The  eigen¬ 
vectors  in  T  will  not,  in  general,  be  orthonormal  nor  even  orthogonal.  The  columns  of  T  are 
orthonormal  if  —  and  only  if  —  H  is  an  orthogonal  matrix  and  orthogonal  if  H  is  a  diagonal 
matrix.  Proof  of  these  statements  will  not  be  made  as  they  amount  merely  to  a  formal  state¬ 
ment  of  the  basic  results  of  that  portion  of  the  theory  of  matrices  which  deals  with  repeated 
eigenvalues  of  a  real  symmetric  matrix. 

A  connection  may  now  be  made  between  the  method  of  this  paper  and  the  textbook  method. 


Assuming  that  the  column  and  coordinate  rearrangements  leading  to  Equation  7  have  been 
carried  out,  one  may  write 


a'  ‘ 

r  a  ! 

B 

a'a 

a'b  • 

E  =  C'  C  - 

B' 

L 

• 

.  B'A 

B'B 

It  follows  that 


a'  a  I  a'  a 

_ _ i _ _ 

B  A  !  B' B 


-A'1  B 
I 


(29) 


(30) 


where  in  Equation  30  the  matrix  on  the  right  is  a  P  by  (P-R)  null  matrix.  It  has  been  previously 
noted  that  the  columns  of  /8  are  linearly  independent.  It  is  quite  clear  from  Equation  30  that 
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the  columns  of  /3  are  eigenvectors  of  E  corresponding  to  P-R  eigenvalues  having  the  value 
zero.  Here  R  represents  the  number  of  rows  of  C  ,  and  by  the  hypothesis  made  in  order  to 
apply  the  textbook  method  R  represents  also  the  rank  of  The  Matrix  E  is  clearly  of 
rank  R  and  therefore  possesses  no  more  than  P-R  eigenvalues  with  value  zero. 

Thus,  the  textbook  solution  is  seen  to  be  a  special  case  of  the  general  solution  which 
would  result  if  the  method  of  this  paper  were  applied  after  the  rearrangements  in  C  leading 
to  Equation  7  were  made. 


SECTION  vn 
FIRST  EXAMPLE 


Here  the  method  is  applied  to  derive  the  equations  of  motion  of  a  simple  chain  of  spring- 
mass  elements.  The  main  intent  is  to  illustrate  application  of  the  method.  However,  some 
points  of  general  interest  will  arise. 


The  system  consists  of  five  point  masses  connected  by  linear  massless  springs  as  shown 
in  Figure  1. 


12  3  4  5 


Figure  1.  Spring-mass  System 


Each  of  the  masses  and  each  of  the  spring  constants  is  assumed  to  have  unit  magnitude.  The 

rr. 

masses  may  displace  only  in  the  horizontal  direction  and  the  displacement  of  the  n  mass  is 
denoted  by  x^.  A  positive  value  of  xnis  taken  to  mean  displacement  to  the  right  and  a  negative 

value  displacement  to  the  left.  A  horizontal  force  F  ,  positive  to  the  right  and  a  function  of 

th  ^ 

time  only,  acts  upon  the  n11  mass. 
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The  five  displacements  constitute  a  set  of  independent  coordinates  which  determine  the 
configuration  of  the  system  at  any  instant;  and  in  terms  of  these  coordinates,  it  is  easy  to 
write  down  directly  equations  of  motion  of  the  system  in  the  form  of  Equation  1  with 

qn=xn  n  =  I,  2  ,  3 , 4,  5  (31  ) 


K  = 


(32) 


(33) 


(34) 


Thus,  from  a  practical  point  of  view,  the  method  of  this  paper  is  not  needed  for  an  analysis  of 
the  system  since  the  end  result  of  the  method,  equations  of  motion  in  terms  of  independent 
coordinates,  is  readily  obtained  by  inspection.  However,  as  an  object  here  is  to  illustrate  the 
method,  let  the  system  be  viewed  in  a  different  way  as  illustrated  inFigure  2.  There  the  system 
of  Figure  1  is  shown  figuratively  divided  into  four  parts  by  cuts  at  the  three  inner  masses, 
producing  an  eight-mass  system. 


#w\i  *v\ A*  I'W'v#  rwv# 

I  2  3  4  5  6  7  8 

Figure  2.  Cut  System 


The  half  circles  represent  masses  of  one-half-unit  magnitude.  The  displacement  of  the  p^1 
mass  of  this  cut  system  is  denoted  by  y  . 
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It  is  assumed  that  three  equations  of  constraint  are  imposed  on  the  coordinates,  namely 

y2  =  y3  (35) 

y4  =  y5  (36) 

yg  =  y?  (37) 

Thus  the  coordinates  y^  are  not  independent,  and  from  the  simple  geometric  considerations 
involved  it  is  clear  that  under  these  equations  of  constraint  the  systems  of  Figure  1  and 
Figure  2  are  the  same.  In  terms  of  the  coordinates  y^  a  Lagrangian  of  the  system  may  be 
expressed  by  an  equation  similar  to  Equation  5  with 

z  -  y  p  =  I,  2,  3,  4,  5,  6,  7,  8  (38) 


The  work  W  for  the  system  may  be  expressed  by  an  equation  following  the  form  of 
Equation  6  with 

Z  =  |  F.  ,0,  F2  ,0,F31  0,F4  ,F6  |  ( 41 } 


It  is  noted  that  the  form  indicated  by  Equation  41  for  the  matrix  Z  is  not  unique.  The  following, 
form,  for  example,  will  serve  equally  well 

Z'  =  |  F,  (l/2)Fz  ,(I/2)F2,(I/3)F3  ,(2/3)F3  ,  (l/5)F4  ,(4/5)F4,FB  (42) 

All  that  is  required  is  that  Z'  when  introduced  in  Equation  6  should  yield  the  work  done  during 
any  displacement  consistent  with  the  equations  of  constraint. 
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Equations  35,  36,  and  37,  the  equations  of  constraint,  may  be  put  in  the  form  of  Equation  4 

with 
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The  first  step  in  the  application  of  the  method  is  to  compute  the  matrix  E  defined  by 
Equation  16.  This  computation  yields 
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-  1 
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0 

0 

0 

<  44 ) 


The  matrix  U  which  follows  is  a  modal  matrix  of  the  Matrix  E  ,  as  may  be  easily  verified 
by  substitution  into  Equations  17  and  18. 


0 

0 

0 

1 

0 

0 

0 

0 

i/yr 

0 

0 

0 

I/-/2- 

0 

0 

0 

-1 A/T 

0 

0 

0 

i/yr 

0 

0 

0 

0 

I//F 

0 

0 

0 

I//T 

0 

0 

0 

-i/yr 

0 

0 

0 

l/7T 

0 

0 

0 

0 

I//2 

0 

0 

0 

1/72 

0 

0 

0 

-i/yr 

0 

0 

0 

1/72 

0 

0 

0 

0 

D 

0 

0 

0 

D 

(4  5) 
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The  Matrix  X  containing  the  eigenvalues  associated  with  the  modal  matrix  is  given  by 


The  first  three  eigenvalues  are  finite  and  the  last  five  have  the  value  zero.  Therefore  the  last 
five  columns  of  U  constitute  a  suitable  transformation  Matrix  T  .  It  follows  that  acceptable 
independent  coordinates  for  describing  any  configuration  of  the  system  consistent  with  the 
equations  of  constraint  are  five  coordinates  qQ  related  to  the  displacements  yp  by  the  equation 


>5 
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*7 

i  y8 
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0 

1 

(47) 


From  the  Equation  Q  =  T'  Z  it  follows  also,  using  either  Equation  41  or  Equation  42,  that 


generalized  forces  suitable  for  use  with  the  coordinates  qn,  are  given  by 


'  0|  ' 

■v 

F, 

°2 

l//2~  F2 

> 

i/ya  f3  - 
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\/Jl  % 

°5 

V.  - 

k  s 

Completing  the  steps  in  the  method  gives 


M  :  T'MT  = 
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Equations  48,  49,  and  50  give  all  the  ingredients  necessary  for  writing  equations  of  motion 

for  the  spring-mass  system  in  the  form  of  Equation  2.  By  use  of  Equation  47  solutions  of  the 

equations  giving  time  histories  of  the  coordinates  qn  can  be  transformed  into  time  histories 

of  the  original  coordinates  yp.  If  initial  conditions  consistent  with  the  equations  of  constraint 

are  given  in  terms  of  the  coordinates  y  the  equations 

P 

,  q  =  T'y  (51) 

and 


q  =  T'y 


may  be  used  to  convert  them  into  initial  conditions  on  the  coordinates  q 


It  may  be  noted  that  the  matrices  K,  M  ,  and  Q  given  in  the  three  preceding  equations  are 
not  identical  to  the  corresponding  matrices  which  were  written  down  directly  from  simple 
physical  considerations  in  Equations  32, 33,  and  34,  respectively.  Either  of  the  sets  of  matrices 
forms  a  valid  basis  for  equations  of  motion  of  the  system  of  Figure  1.  The  difference  between 
the  matrices  arises  from  the  fact  that  the  coordinates  qn  arrived  at  by  the  method  of  this 
paper  are  not  related  to  the  coordinates  yp  in  the  same  way  as  are  the  coordinates  x  . 
Equation  47  shows  the  relationship  between  the  coordinates  qn  and  y  whereas  coordinates 
xn  yp  are  related  by  the  equation  ^ 
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Equation  53  may  be  written 


where 


y  =  Tx 

T  =  TH 


(54) 

(55) 


in  which  T  is 
and 


the  transformation  matrix  in  Equation  47 ,  derived  by 
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JT 
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1 

the  method  of  this  paper, 


(56) 


Thus  the  coordinates  x^,  which  represent  displacements  of  masses,  are  in  the  category  of 
coordinates  discussed  in  connection  with  Equation  27.  The  foregoing  discussion  illustrates  a 
feature  of  the  method  of  this  paper  which  should  be  recognized  by  anyone  using  the  method. 
That  is:  The  coordinates  qfl  produced  by  the  method  are  generally  abstract  in  character  and 
do  not  lend  themselves  to  simple  physical  interpretations. 

It  is  instructive  to  reexamine  the  matrix  C  in  Equation  43  and  think  about  the  decisions 
involved  in  applying  the  textbook  method.  Consider  the  three  pairs  of  displacements  (y2«y3>. 
(y  ,y  ),  and  (y6,y7)  straddling  the  cuts  in  Figure  2.  Let  triplets  of  displacements  be  formed 
by  taking  one  and  only  one  displacement  from  each  pair,  for  example  (y2,y5»y7)  and  (y3»y4»y6)- 
If  the  displacements  in  any  such  triplet  are  taken  to  make  up  the  elements  of  the  column  z*a) 
in  Equation  7  the  Matrix  A  formed  from  the  corresponding  columns  will  be  nonsingular  and 
the  textbook  method  will  succeed.  If  the  elements  of  z  are  chosen  from  among  the  eight 
coordinates  y  in  any  otherway,  the  Matrix  A  will  be  singular.  In  applying  the  textbook  method 
to  this  simple  problem,  recognition  of  the  combinations  of  coordinates  suitable  to  form  z 
must  come  about  either  from  physical  insight  or  from  understanding  of  linear  dependence 
among  the  columns  of  C  .  In  applying  the  method  of  this  paper  it  is  not  necessary  to  think 
directly  about  the  physics  or  about  linear  dependence.  Instead,  the  problem  becomes  one  of 
finding  a  modal  matrix  of  E  and  identifying  the  columns  associated  with  eigenvalues  having 
the  value  zero.  Due  to  the  block-diagonal  form  of  E  in  this  case,  it  was  possible  by  inspection 
to  put  down  exactly  a  modal  matrix  and  the  eigenvalues  of  E  .  Therefore,  all  decisions  in  ap¬ 
plication  of  the  method  of  this  paper  could  be  made  easily  on  a  purely  mathematical  basis. 
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SECTION  vm 
SECOND  EXAMPLE 


The  purpose  here  is  to  discuss  an  example  in  which  redundancies  in  the  equations  of 
constraint  arise  in  a  natural  way.  The  mechanical  system  is  shown  in  Figure  3.  A  cylindrical 
elastic  shell  is  fixed  at  one  end 


to  an  immovable  base.  At  the  other  end  a  thin  massive  rigid  disk  is  attached  to  the  wall  of 
the  shell  by  four  symmetrically  placed  pins.  Points  in  the  shell  wall  are  assumed  to  displace 
only  longitudinally. 

Adopting  an  approximation  common  in  practical  vibration  analysis,  the  longitudinal 
displacement  u  of  a  general  point  in  the  shell  wall  is  expressed  as  a  linear  combination  of  a 
finite  number  of  displacement  functions.  The  expansion  assumed  is 

“  :  I  I  8  +8  cos  40]  sin  m  (57) 

m,o  m  ,4  2  i 

where  m  takes  on  integral  values  and  the  summation  sign  indicates  summation  of  the  terms 

corresponding  to  some  finite  number  of  selected  values  of  m.  The  coefficients  8  and 

a  •  m*° 

®  m>4  are  functions  of  time  alone  and  serve  as  coordinates  which  describe  the  instantaneous 

configuration  of  the  shell. 


325 


AFFDL-TR-68-150 


Assuming  small  displacements,  the  instantaneous  position  of  the  disk  is  determined  by 
specification  of  three  coordinates  8^,  a  £  ,  and  a  ^  defined  as  follows: 


(1)  8  is  the  displacement  of  the  center  of  the  disk  parallel  to  the  longitudianal  axis  of 


the  shell. 


(2)  ct£  and  a  ^  are  small  rotations  about  axis  £  and  17  ,  respectively,  as  shown  in  the 
sketch. 


Equating  the  displacements  of  the  disk  to  the  displacements  of  the  shell  at  each  of  the 
four  pins  gives 

8C  -RC£  =  I  (Sm  o  +  Sm  4)sin 


8C  +Ra  =  I(8m  +  8m  .)sin^ 


m.o  m,4 


8  +  Raf  =  X  (S_  +8  )  sin 

c  f  m,o  m,  4  d 


mw 


8.  —  Ra_  =  X  (S_  „  +  8_  . )  sin  _ 
c  7j  m,o  m,  4  2 


(58) 


where  R  is  the  radius  of  the  cylinder.  If,  in  the  summation  on  the  right,  only  the  terms 
corresponding  to  m  -  1  are  retained,  the  equations  may  be  put  in  the  form 


where 


and 


CZ:0 


(59) 


(60) 


(6!) 


It  will  be  clear  on  inspection  that  an  attempt  to  arrive  at  independent  coordinates  for  this 
system  by  a  straightforward  application  of  the  textbook  method  must  fail  because  any  choice 
of  the  Matrix  A  will  lead  to  a  matrix  which  has  at  least  two  identical  columns  and  which  is 
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therefore  singular.  This  difficulty  stems  from  the  fact  that  the  system  of  equations  is  re¬ 
dundant  which  may  be  demonstrated  by  adding  Rows  one  and  three  of  Matrix  C  and  subtracting 
from  the  result  Row  two  producing  Row  four. 

One  way  to  arrive  at  independent  coordinates  would  be  to  discard  the  fourth  equation  from 
the  system  and  apply  the  textbook  method  to  the  first  three  equations.  However,  this  approach 
requires,  in  general,  the  following: 

(1)  Recognition  in  the  first  place  that  the  system  is  redundant. 

(2)  Recognition  of  dependent  equations. 

(3)  Recognition  of  a  nonsingular  submatrix  A  after  redundant  equations  are  discarded. 

For  the  example  problem  under  consideration,  the  required  understanding  of  the  structure 
of  the  equations  may  be  gained  by  inspection.  In  practical  work,  however,  there  may  be  many 
equations  of  constraint  involving  many  unknowns,  and  the  coefficients  making  up  the  Matrix  C 
will  usually  not  be  small  integers.  Generally,  in  such  situations,  little  of  use  can  be  deduced 
about  the  system  merely  by  inspection  of  the  matrix  of  coefficients.  Also,  one  cannot  always 
rely  on  physical  insight  to  detect  and  understand  redundancies.  In  fact,  in  the  example  being 
considered,  there  is  nothing  on  the  face  of  it  in  the  physics  to  warn  of  a  redundancy.  Further, 
there  are  considerable  theoretical  and  practical  difficulties  in  making  computational  tests  for 
singularity  and  redundancy  when  there  is  error,  such  as  roundoff  error,  in  the  process  by 
which  the  coefficients  of  the  equations  of  constraint  are  generated. 


Proceeding  now  to  apply  the  method  of  this  paper,  the  Matrix  E  is  given  by 
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The  eigenvalues  of  E  are 

X|=I21X2:2,X3  =  2)X4=0,\b=0  (62) 
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It  may  be  easily  verified  that  the  two  columns  of  the  Matrix  T  which  follows  are  orthonormal 
eigenvectors  of  E  corresponding  to  the  two  eigenvalues  ^  and  ^  ^  which  have  the  value  zero. 


l//6“ 

i/yr 

0 

0 

0 

0 

-i/ye" 

+  i/y“T 

+  2//6“ 

0 

(63) 


Therefore  the  system  may  be  described  by  two  independent  coordinates  and  q2  related  to 
the  coordinates  in  z  by  the  equation 


z  --  T  q 


(64) 


As  can  be  seen,  direct  concern  with  the  number  and  nature  of  redundancies  in  the  equations 
of  constraint  is  unnecessary  when  the  method  of  this  paper  is  used.  The  problem  reduces  in 
substance  to  that  of  determining  an  orthonormal  set  of  eigenvectors  of  E  corresponding  to 
the  eigenvalues  of  E  which  have  the  value  of  zero. 
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SECTION  IX 

COMMENTS  ON  NUMERICS 


In  the  examples  it  was  possible  to  put  down  exactly  the  Matrix  c  ,  to  carry  out  exactly  the 
the  multiplication  C'C  producing  the  Matrix  E  ,  and  to  determine  exactly  the  eigenvalues  of  E 
and  orthonormal  eigenvectors  corresponding  to  the  eigenvalues  with  value  zero.  In  practical 
work,  however,  numerical  error  due  to  roundoff  and/or  truncation  may  be  introduced  at  any  of 
these  three  stages  of  calculation.  The  extreme  effect  of  such  errors  would  be  complete  loss 
of  numerical  significance  in  the  digits  representing  the  eigenvalues  of  E  and  the  elements  of 
the  eigenvectors  of  E .  In  the  event  the  computation  is  subject  to  serious  loss  of  significance, 
the  Matrix  C  is  said  to  be  ill-conditioned  with  respect  to  the  computing  process  used.  The 
best  indication  of  ill-conditioning  is  sensitivity  of  final  results  to  small  changes  in  the  elements 
of  C .  The  authors  have  applied  the  method  of  this  paper  a  number  of  times  in  practical 
vibration  analysis  and  have  not  as  yet  encountered  a  situation  in  which  the  Matrix  C  is 
ill-conditioned.  Speaking  from  general  experience,  however,  the  possibility  of  ill-conditioning 
must  be  anticipated  whenever  simultaneous  equations  are  solved  numerically,  and  the  method 
of  this  paper  presents  no  exception  to  this  statement.  When  an  ill-conditioned  system  arises, 
the  recourse  most  often  open  is  to  increase  the  number  of  digits  carried  in  the  computation. 
If  this  is  attempted  in  connection  with  the  method  of  this  paper,  it  should  be  recognized  that 
it  may  be  necessary  to  increase  the  carried  significance  in  the  stage  of  the  calculation  in 
which  the  elements  of  C  are  generated  as  well  as  in  the  implementation  of  the  multiplication 
c'  C  and  in  calculating  the  eigenvalues  and  eigenvectors  of  E. 

Another  consequence  of  numerical  error  is  that  finite  numbers  may  be  generated  for 
eigenvalues  of  E  which  would  be  precisely  zero  if  there  were  no  error  in  the  computing 
process.  This  raises  the  question,  in  principle  at  least,  of  the  possibility  of  rigorous  dis¬ 
tinction  between  finite  numbers  representing  finite  eigenvalues  of  E  and  finite  numbers 
representing  eigenvalues  of  E  which  are,  in  fact,  zero.  In  the  authors’  experience  this  has 
proved  to  be  more  of  a  problem  in  principle  than  in  practice.  The  authors  use  the  threshold 
Jacobi  Method  to  compute  the  eigenvalues  and  a  modal  matrix  of  E .  Approximately  15 
significant  figures  are  carried  throughout  the  calculation.  With  this  procedure,  inspection  of 
the  eigenvalues  computed  for  E  has  always  revealed  two  clearly  distinguishable  sets  of 
numbers,  the  numbers  in  one  set  being  many  orders  of  magnitude  smaller  than  the  numbers 
in  the  other.  The  set  of  numbers  with  relatively  large  magnitudes  are  regarded  as  finite 
eigenvalues,  and  the  remaining  numbers  are  considered  to  be  eigenvalues  with  value  zero. 
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It  is  helpful  to  recognize  that  the  number  of  finite  eigenvalues  of  E  can  be  anticipated 
if  the  rank  R  of  C  is  known.  For  the  rank  of  E  is  equal  to  R,  and  it  is  not  difficult  to  show  that 
the  number  of  finite  eigenvalues  of  E  is  therefore  equal  to  R.  Frequently,  it  is  known  from 
physical  or  geometric  considerations  that  the  equations  of  constraint  are  linearly  independent 
in  which  case  the  rank  R  of  C  is  equal  to  the  number  of  rows  of  C . 


As  has  been  indicated,  in  the  authors’  experience  with  the  method  under  discussion,  it 
has  always  been  possible  to  distinguish  with  confidence,  on  the  basis  of  magnitude,  between 
numbers  representing  finite  eigenvalues  and  finite  numbers  representing  eigenvalues  which 
are,  in  fact,  zero.  If  such  a  distinction  could  not  be  made,  that  is,  if  the  eigenvalues  of  E 
were  to  decrease  gradually  from  maximum  to  minimum,  the  authors  would  suspect  ill- 
conditioning  of  the  Matrix  C  . 


SECTION  X 

CONCLUDING  REMARKS 

A  computational  method  has  been  devised  by  which  equations  of  motion  of  a  linear  me¬ 
chanical  system  in  terms  of  independent  coordinates  can  be  generated  when  basic  information 
about  the  system  is  available  in  terms  of  coordinates  which  are  not  independent  but,  instead, 
are  governed  by  linear  homogeneous  equations  of  constraint.  Necessity  for  this  derivation 
arises  frequently  in  practical  vibration  analysis.  The  method  is  believed  to  be  new,  and  ex¬ 
perience  in  analyzing  the  vibrations  of  shells  indicates  that  it  will  very  often  offer  decided 
advantages  over  methods  previously  used.  In  the  method,  a  real  symmetric  matrix  is  con¬ 
structed  by  an  operation  which  involves  only  the  coefficients  in  the  equations  of  constraint. 
The  eigenvectors  and  corresponding  eigenvalues  of  the  symmetric  matrix  are  computed. 
Then,  a  transformation  matrix  leading  to  independent  coordinates  is  assembled  from  the 
eigenvectors  corresponding  to  eigenvalues  having  the  value  zero.  The  main  advantages  of  the 
method  are:  (1)  For  the  most  general  constraint  equations,  the  problem  is  reduced  to  calcu¬ 
lating  the  eigenvectors  and  eigenvalues  of  a  symmetric  matrix.  This  calculation  is  one  of  the 
most  successful  applications  of  modern  digital  computers;  (2)  The  method  is  applicable  when 
there  are  redundancies  in  the  equations  of  constraint. 
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